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We model the initial confinement-deconfinement transition in relativistic heavy-ion collisions as 
a rapid quench in view of expected rapid thermalization to a QGP state. The transition is studied 
using the Polyakov loop model, with the initial field configuration (in the confining phase) covering 
a small neighborhood of the confining vacuum I ~ 0, as appropriate for T < T c . Quench is imple- 
mented by evolving this initial configuration with the effective potential at a temperature T > T c . 
We study the formation of Z(3) domain structure and its evolution during the transition as I rolls 
down in different directions from the top of the central hill in the effective potential of I. When 
explicit Z(3) symmetry breaking effects (arising from dynamical quark effects) are small, then we 
find well defined Z(3) domains which coarsen in time. Remarkably, the magnitude plot of I shows 
vacuum bubble like configurations arising during the quench. This first order transition like behav- 
ior occurs even though there is no metastable vacuum separated by a barrier from the true vacuum 
for the parameter values used. When the initial field configuration everywhere rolls down roughly 
along the same direction (as will happen with large explicit symmetry breaking) then we do not 
find such bubble-like configurations. However, in this case we find huge oscillations of I with large 
length scales. We show that such large oscillations can lead to large fluctuations in the evolution of 
flow anisotropies compared to the equilibrium transition case. 
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I. INTRODUCTION 



In relativistic heavy-ion collision experiments (RHICE) the collision of two nuclei leads to a hot dense region which 
is expected to rapidly achieve state of thermal equilibrium. For the relevant range of energies and colliding nuclei 
at RHIC and at LHC there is compelling evidence that a region of quark-gluon plasma is created in these collisions. 
Simulation results as well as experimental data, such as elliptic flow measurements, all point towards a very rapid 
thermalization to the QGP phase, within a proper time less than 1 fm. We thus have a system which starts out 
in the confining phase, and within proper time of (probably much less than) 1 fm makes a transition to the QGP 
phase (with maximum temperature estimates ranging from 200 MeV to more than 700 MeV for the relevant energies 
in these experiments). Lattice results show that for real world QCD with small baryon density (as appropriate for 
the central rapidity regions in RHICE) the transition is likely to be a crossover. With that, the dynamics of the 
transition depends crucially on the rate of temperature change compared to the time scale of the evolution of the 
order parameter field. For an equilibrium transition, we had studied the formation of Z(3) domains and associated 
QGP strings using a first order transition dynamics via bubble nucleation [312]. The transition was simulated using 
the Polyakov loop, l(x), as an order parameter for the confinement-deconfinement transition [3J. These studies were 
appropriate for large chemical potential case, as in lower energy collisions, where the transition is expected to be of 
first order. Though, the results for Z(3) wall network etc. having certain universal characters, may be applicable in 
more general context, even for a cross-over, as explained in [IJH]. 

However, given the very short time scale of initial thermalization to QGP state, an equilibrium dynamics of the 
transition appears unlikely. A more appropriate description of the transition should employ quenched dynamics in 
which the growth of Z(3) domains will be via spinodal decomposition. In this paper we carry out such a study using 
the Polyakov loop, l(x), as an an order parameter for the confinement-deconfinement transition, with an effective 
potential of the kind used in refs.[3]-[5]. For our simulation results we choose a definite parametrization of the effective 
potential [5] as was used in our previous works [IJ [2]. We model the phase transition in this Polyakov loop model, 
with the initial field configuration (in the confining phase) covering a small neighborhood of the confining vacuum 
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I ~ as appropriate for the initial T = system. In a quench, the temperature rapidly (rather suddenly) increases to 
its maximum value To with the effective potential changed accordingly. The initial field I, unable to relax to the new 
equilibrium vacuum state in this short time, becomes unstable and rolls down in different directions from the top of 
the central hill in the effective potential of I. We study the formation of Z(3) domain structure during this evolution. 
When explicit Z(3) symmetry breaking effects (arising from dynamical quark effects) are small, then we find well 
defined Z(3) domains which coarsen in time. With a symmetric initial patch of I, all the three Z(3) domains form 
with random shapes and rapidly increase in size by coarsening. Remarkably, the magnitude plot of I shows vacuum 
bubble like configurations, such as those which arise in a first order transition, arising during the quench in this case 
(when the initial field rolls down in different directions). This first order transition like behavior occurs even though 
there is no metastable vacuum separated by a barrier from the true vacuum for the parameter values used. These 
bubble like configurations expand as well, somewhat in similar manner as during a first order transition. 

When the initial patch of I is only partially symmetric around I = (as appropriate for small explicit symmetry 
breaking from quark effects), the dynamics retains these qualitative aspects with expected changes. Thus, true vacuum 
domains (with 9 — 0) are more abundant and they also grow fast at the cost of the other two Z(3) domains (which are 
now metastable due to explicit symmetry breaking) . Still, for small explicit symmetry breaking, all the three types of 
domains occur during the initial stages of the transition. There are few, small Z(3) walls as there are fewer metastable 
Z(3) domains embedded in the dominant true vacuum with 9 — 0. These walls shrink rapidly and eventually only the 
true vacuum survives. 

The dynamics is found to be very different when the explicit symmetry breaking due to quark effects is taken 
to be strong. In this case the initial patch of I (around equilibrium point for T = effective potential) could be 
significantly shifted towards the true vacuum for the quenched T = Tq effective potential. In such a situation, I will 
roll down roughly along the same direction with angular variations becoming smaller during the roll down. In this 
case only 9 = vacuum survives and no other Z(3) domains are formed. Also, in this case we do not find bubble-like 
configurations. However, in this case we find huge oscillations of I with large length scales. This behavior is known 
from previous studies of the dynamics of scalar field in a quench [7] and is expected here when angular variations 
are small. The dynamics of field in such a case is dominated by large length scale coherent oscillations. It leads to 
novel scenarios of reheating via parametric resonance in the case of inflation in the early Universe [? ] . In our case of 
RHICE also it raises important questions about the possibility of parametric resonance and of novel modes of particle 
production from these large oscillations of I during the early stages of the transition. In the present work we explore 
another important effect of these large oscillations, on the evolution of flow anisotropies in RHICE. As most of the 
flow anisotropies are expected to develop during first few fm of the QGP formation, it becomes an important question 
if the presence of large oscillations of I can affect the development of these flow anisotropies. As we will see, this 
indeed happens and these large oscillations lead to large fluctuations in flow anisotropy. 

We mention here that we do not discuss which of the cases discussed here (small, or large explicit symmetry 
breaking) may actually be realized in RHICE. This is primarily because of lack of understanding of the magnitude of 
explicit symmetry breaking term near T = T c . (Some discussion of this has been provided in our earlier works [lj[2].) 
Also, the spread of initial field configuration about confining vacuum plays a crucial role here and that in turn depends 
on details of pre-equilibration stage. A proper understanding (or modeling) of this stage and resulting estimate of the 
initial spread is essential before one can make more definitive statements about the dynamics of the transition. We 
hope to discuss these issues in future. Our main purpose here is to illustrate the possibility of very different types 
of dynamics of transition depending on the initial configuration. Our results show that in quenched dynamics, Z(3) 
domains, and resulting Z(3) domain walls, can last for any reasonably length of time only when explicit symmetry 
breaking terms are very small. Otherwise either different domains do not form at all, hence no Z(3) walls are formed, 
or fewer metastable Z(3) domains form embedded in most abundant true vacuum region. In the latter case resulting 
Z(3) walls are smaller to begin with, and disappear quickly. 

Importance of Z(3) walls in RHICE has been discussed by us in previous works, where we have also emphasized 
non-trivial scattering of quarks from Z(3) walls. We have explored its consequences for cosmology as well as for 
RHICE [U El HI [9] , including the possibility of CP violating scattering of quarks from Z(3) walls leading to interesting 
observational implications [10] - Recently an interesting possibility has been discussed by Asakawa et al. in ref. 
[llj where it is argued that scattering of partons from Z(3) walls may account for small viscosity as well as large 
opacity of QGP. Our results for the formation of Z(3) domains during quench (which are abundant only for small 
explicit symmetry breaking case) can be relevant for the studies in ref. [llj . In this context we mention that smallest 
reasonable size Z(3) domains (hence Z(3) walls) we find are of order 1-2 fm at very early times, and at that stage the 
magnitude of the Polyakov loop order parameter I is very small, of order few percent of its vacuum expectation value. 
The quark scattering from Z(3) walls is likely to be small for such a small magnitude of / [51-HU]. By the time the 
magnitude of I becomes significant, domains coarsen to have large sizes, of order several fm. Thus, in the context of 
our model it appears difficult to form very small Z(3) domains which still can scatter partons effectively (as needed 
in the study of ref.[TT]). 
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It is important to note that our results for domain growth and fluctuations for the quench case are dominated by 
the spinodal instabilities during the roll down of the field from the top of the hill of the effective potential. These 
instabilities arise primarily from the nature of the quench when the initial field configuration becomes unstable to 
exponential growth of long wavelength modes due to sudden change in the shape of the effective potential (the quench) , 
and will be in general present even if the transition is a cross-over. In this sense we believe that the qualitative aspects 
of our results have a wider applicability and are not crucially dependent on the specific form of the effective potential 
used here. 

We briefly mention here that issues related to the physical meaning of Z(3) domains etc. have been discussed in 
the literature (see, discussion of these in our earlier works [TJ [2] and we will not repeat it here. We only note that 
recent work of Deka et al. [12] has provided a support for the existence of these metastable vacua from Lattice. 
We also note that a simulation of spinodal decomposition in Polyakov loop model has been carried out in ref. [13] , 
where fluctuations in the Polyakov loop and growth of long wavelength modes (representing domain formation) are 
investigated. In comparison, the main focus of our work is on detailed growth of domains due to coarsening with 
and without explicit symmetry breaking, new bubble like structures, and existence of large fluctuations affecting flow 
anisotropies in important ways. 

The paper is organized in the following manner. In section II, we discuss the essential aspects of the Polyakov loop 
model of confinement-deconfinement phase transition and the effective potential used from ref. ;5 . . Section III presents 
the numerical technique of simulating the phase transition via quench. We first discuss the case without any explicit 
symmetry breaking and study the formation and evolution of Z(3) domains during the quench. Sec. IV discusses 
the unexpected result of bubble like structures arising during the quench. We also discuss the case of small explicit 
symmetry breaking term when true vacuum domains become more abundant, but overall picture of the transition 
remains roughly similar. In Sec.V we discuss the case when explicit symmetry breaking effects are strong leading to 
I rolling down everywhere roughly along 9 — 0. This leads to large oscillations of I. In Sec. VI we study the effects 
of these large oscillations on flow anisotropy and show that it leads to large fluctuations in elliptic flow and spatial 
eccentricity (as compared to the case of equilibrium transition). Section VIII presents conclusions. 



II. MODELING THE PHASE TRANSITION 



We briefly recall the salient features of the model we use for the confinement-deconfinement phase transition. The 
order parameter is taken to be the the expectation value of Polyakov loop l(x), 

1 

l(x) = -Tr{Vexp{ig J A (x,r)dr)) (1) 

Where A^(x : t) is the time component of the vector potential A^(x,t) = A®(x,T)T a , T a are the generators of 
SU(N) in the fundamental representation, V denotes path ordering in the Euclidean time r, g is the gauge coupling, 
and (3 — l/T with T being the temperature. N (— 3 for QCD) is the number of colors. The complex scalar field l(x) 
transform under the global Z(N) (center) symmetry transformation as 



l(x) -» exp{2nin/N)l{x), n = 0, 1, ...(N - 1) (2) 

The expectation value of l(x) is related to e~$ F where F is the free energy of an infinitely heavy test quark. For 
temperatures below T c , in the confined phase, the expectation value of Polyakov loop is zero corresponding to the 
infinite free energy of an isolated test quark. (Hereafter, we will use the same notation l(x) to denote the expectation 
value of the Polyakov loop.) Hence the Z(N) symmetry is restored below T c . Z(N) symmetry is broken spontaneously 
above T c where l(x) is non-zero corresponding to the finite free energy of the test quark. Effective theory of the 
Polyakov loop has been proposed by several authors with various parameters fitted to reproduce lattice results for 
pure QCD [3H5]. We use the Polyakov loop effective theory proposed by Pisarski [5,6. The effective Lagrangian 
density can be written as 

L=^\d,l\ 2 T 2 -V(l) (3) 
Where the effective potential V(l) for the Polyakov loop, in case of pure gauge theory is given as 

V(l) = (^\l\ 2 - b j(l 3 + (D») + \{\l\ 2 ) 2 W (4) 
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At low temperature where 1=0, the potential has only one minimum. As temperature becomes higher than T c the 
Polyakov loop develops a non vanishing vacuum expectation value Iq, and the cos39 term, coming from the I 3 + l* 3 
term above leads to Z(3) generate vacua. Now in the deconfined phase, for a small range of temperature above T c , the 
/ = cxtremum becomes the local minimum (false vacuum) and a potential barrier exist between the local minimum 
and global minimum (true vacuum) of the potential. (However, the quench temperature used here is much higher 
than this range and there is no barrier present at such temperatures.) 

The effects of dynamical quarks is included in terms of explicit breaking of the Z(3) symmetry which is represented 
in the effective potential by inclusion of a linear term in I [5J [T3] . The potential of Eq. (4) with the linear term 
becomes, 

v(i) = + n - |i*i 2 - |a 3 + r 3 ) + \m 2 ) 2 )b^ ( 5 ) 

Here coefficient b\ measures the strength of explicit symmetry breaking. The coefficients b±, &2> ^3 and 64 are 
dimensionless quantities. With 61 = 0, other parameters 62, ^3 and 64 are fitted in ref.0 E] such that that the 
effective potential reproduces the thermodynamics of pure SU(3) gauge theory on lattice [T4HI6] . The values of 
various coefficients from refs. [5J |B] are the same as used in our previous works [TJ [5] (including discussions about 
explicit symmetry breaking strength bi) and we do not repeat that discussion here. (With those values of parameters, 
the transition temperature is taken to be T c = 182 MeV.) 



III. NUMERICAL SIMULATION 



In this work, we carry out a 2+1 dimensional field theoretic simulation of the dynamics of confincmcnt-dcconfinement 
transition in a quench. First we work within the framework of Bjorken's boost invariant longitudinal expansion model 
|17j (without any transverse expansion) for the central rapidity region in RHICE. To model the quench, we take the 
initial field configuration to constitute a small patch around I = which corresponds to the vacuum configuration 
for the confining phase. This is for the case of zero explicit symmetry breaking. We will discuss the case of explicit 
symmetry breaking later. We have taken the initial phase of I to vary randomly between and 2ir from one lattice 
site to the other, while the magnitude of I is taken to vary uniformly between and e. (We have also taken the initial 
magnitude to have a fixed value equal to e and the results are similar.) Value of e is taken to be much smaller than 
the vacuum expectation value (vev) of I at T = Tq and results remain qualitatively the same for similar small values 
of e. We report results for e = 0.01 times the vev of I. We take the quench temperature T = 400 MeV. 

This initial field configuration, which represents the equilibrium field configuration of a system with T < T c , is 
evolved using the effective potential with T = Tq > T c . This represents the transition dynamics of a quench. The 
field configuration is evolved by time dependent equation of motion in the Minkowski space [18] 



dr 2 t 8t dx 2 dy 2 2NT 2 dh ' J { ' 



with -g^ = at t = and I = l± + ifa. The temperature is taken to decrease starting with the value To (at an 
initial time r which we take to be 1 fm in all the simulations) as r -1 / 3 as appropriate for the longitudinal expansion 
model. We mention that later in Sect.V we will consider an isotropic geometry for the transverse dynamics of QGP 
as relevant for RHICE. There we will also model non-zero transverse expansion and then the central temperature will 
decrease faster than t -1 / 3 . Here, we take a 2000 x 2000 square lattice with physical size 20 fm x 20 fm. We take this 
lattice as representing the transverse plane of the QGP formed in a central collision and consider mid rapidity region. 
The evolution of field was numerically implemented by a stabilized leapfrog algorithm of second order accuracy both 
in space and in time with the second order derivatives of U approximated by a diamond shaped grid [2] . We evolve 
the field using the periodic, fixed, and free boundary conditions for the square lattice. Here we present the results 
with periodic boundary conditions. We use Ax = 0.01 fm for the present case, later on we use different values of 
Ax, as shown by the lattice size in the corresponding figures, we take At = Ax/y/2 as well as At = 0.9Ax/V2 to 
satisfy the Courant stability criteria. The stability and accuracy of the simulation is checked using the conservation 
of energy during simulation. The total energy fluctuations remains few percent without any net increase or decrease 
of total energy in the absence of dissipative I term in Eq.(6). 

First, we present results for the formation and evolution of Z(3) domains. The initial field configuration in the 
neighborhood of I — becomes unstable when evolved with the effective potential (Eq.(5)) with T = Tq = 400 MeV 
(at t = t = 1 fm). As / rolls down in different direction, settling in one of the three Z(3) vacua, different Z(3) 
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domains form. Initially, as the phase of I is taken to vary randomly from one lattice site to the next, there are no well 
defined domains. Also, the magnitude of I is very small initially making any association of Z(3) structure meaningless 
at such early stages. Situation remains similar for very early times as seen in Fig. la at an early stage r = 1.2 fm (i.e. 
0.2 fm after the quench). In Fig.l we have shown the values of the phase of I around the three Z(3) vacua in terms 
of different shades (colors) to focus on the evolution of Z(3) domain structure. Thus all the values of the phase 9 of I 
are separated in three ranges, between — 2tt/6 to 2ir/6 (9 = vacuum), between 2-7r/6 to tt (9 = 2ir/3 vacuum), and 
between it to 2ir — 2ir/6 (9 = 4ir/3 vacuum). As the field magnitude grows, the angular variation of I also becomes less 
random over small length scales, leading to a sort of Z(3) domain structure. Z(3) domains become more well defined, 
and grow in size by coarsening as shown in sequence of figures Fig.lb-d. Different shades (colors) in Fig.l represent 
the three Z(3) domains. Fig.lb-d shows the growth of these domains at r = 2.0, 2.4, and 2.8 fm, as I relaxes to the 
three Z(3) vacua and domains grow in size by coarsening. The magnitude of I is about 0.04, 0.08, 0.2, and 0.4 for 
Figs.la,b,c, and d respectively. Note that domains grow rapidly to size of order 2 fm within a time duration of about 
1 fm as shown in Fig. lb. Within another 1 fm time, domain size is about 4 fm as seen in Fig. Id. 

The boundaries of different Z(3) domains represent Z(3) walls, and junction of three different Z(3) domains gives 
rise to the QGP strings. These objects have been discussed in detail in our earlier works [Tl [21 119j. 





FIG. 1: 

(a) Evolution of Z(3) domain structure during a quench. Note we take the initial temperature to be To (= 400 MeV) at proper 
time r = ro = 1 fm to set the initial conditions for the simulation. Different shades (colors) represent the three Z(3) domains, 
(a) - (d) show the growth of these domains at t = 1.2, 2.0, 2.4, and 2.8 fm, (with corresponding values of temperature T = 
376, 317, 298, 283 MeV) as I relaxes to the three Z(3) vacua and domains grow in size by coarsening. The magnitude of I is 
about 0.04, 0.08, 0.2, and 0.4 for (a), (b), (c), and (d) respectively. Note that domains grow rapidly to size of order 2 fm 
within a time duration of about 1 fm as shown in (b). Within another 1 fm time, domain size is about 4 fm (as seen in (d)). 

We mention here that smallest reasonable size Z(3) domains (hence Z(3) walls), which we find in our simulation, 
are of order 1-2 fm at very early times, such as seen in Fig. lb. At this stage, the magnitude of the Polyakov loop 
order parameter I is still very small, of order few percent of its vacuum expectation value. This is important when 
one considers the possibility of nontrivial scattering of partons from Z(3) walls [8Ul 1J . The quark scattering from Z(3) 
walls is likely to be small for such a small magnitude of I [8TU0j. By the time the magnitude of I becomes significant, 
domains coarsen to have large sizes, of order several fm, as in Fig.c,d where the magnitude of I is about 20% and 
40%, respectively, of its vacuum expectation value. Thus, in the context of our model, with quenched dynamics of 
transition, it appears difficult to form very small Z(3) domains which still can scatter partons effectively (as needed 
in the study of ref.|llj). 
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IV. BUBBLE LIKE STRUCTURES DURING QUENCH 



In this section we show a very unexpected result. Fig. 2 shows sequence of plots of the magnitude of I during the 
quench described in the previous section. We note the appearance of bubble like structures. These structures are 
also seen to grow in a manner similar to the bubbles for a conventional first order phase transition, as in ref. [TJ [5] . 
However, the distribution of the phase of I does not show any specific local variation related to these bubble-like 
configurations. In a roughly uniform region of the phase these localized bubble-like configurations arise and expand. 
We show here the plots of / in Fig. 2 corresponding to the initial central temperature T = T = 500 MeV (at r = r = 
1 fm). The plots in Fig.2a,b,c,d are for r = 2.0, 2.8, 3.0, and 3.3 fm respectively, with the corresponding values of the 
central temperature being T — 395, 355, 345, 336 MeV. It is important to note that at such high temperatures there 
is no metastable confining vacuum at I = in the effective potential [2]. (Metastable confining vacuum exists from 
T = T C = 182 MeV upto T ~ 250 MeV.) There is no tunneling modeled, (or thermal hopping over the barrier) here, 
nor is expected. One expects a simple roll down of the field representing the dynamics of spinodal decomposition 
during the quench. We mention that similar bubble-like configurations also arise with To = 400 MeV. However, in 
that case the temperature range goes below T — 250 MeV. To make a clear case that these bubble like configuration 
have nothing to do with a first order transition like situation (even remotely), we have shown plots with T = 500 



These plots correspond to initial central temperature To = 500 MeV. (a)-(d) show the surface plots of |/| during early stages 
with subsequent formation and growth of bubble like structures just like a first order phase transition. Plots in (a) -(d) are at 
t — 2, 2.78, 3 and 3.3 fm respectively, with corresponding temperature being T — 395, 355, 345, 336 MeV. These bubble-like 
configurations are surprising as there is no barrier here, and no metastable vacuum in the effective potential for this 

temperature range. 

This result is very unexpected and points to new interesting possibilities for the phase transition dynamics. For 
example, such bubble like structures may lead to a dynamics of phase separation in the case of the universe similar to 
the original Witten's scenario [20] . even when there is no underlying first order transition. This may also have impor- 
tant implications for RHICE. More importantly this new possibility of transition dynamics needs to be understood 
and analyzed in detail, see [3T] for a study of these issues. 

The studies of this section and the previous section apply to the case without any explicit symmetry breaking, i.e. 
with bi = in Eq.(5). We have repeated these simulations with small explicit symmetry breaking effects on the initial 
conditions. We take b\ = 0.005 as in ref. [2] here as well as in the next section. By small explicit symmetry breaking 
effects on the initial conditions we mean that the initial patch of I is taken to shift towards 9 = vacuum for T = Tq 
effective potential, while still overlapping with the initial equilibrium value of I. This forces I to roll down to different 
9 directions at least at some fraction of lattice points, though major fraction now rolls down towards 9 = 0. Fig. 3 
shows sequence of plots showing growth of Z(3) domains in one such case. We see that one of the vacua (9 = 0) 



MeV. 





FIG. 2: 
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expands dominantly while other domains remain relatively smaller. The Z(3) domain walls in this case are smaller 
and disappear faster compared to the case without explicit symmetry breaking. The dynamics of bubble like structure 
retains its qualitative aspects in this case as long as the field rolls down in different directions. 




10 -8 -6 -4 -2 2 4 6 8 10 -10 -8 -6 -4 -2 2 4 6 8 10 

?c (fin) ?i (fin) 



FIG. 3: 

(a) Field configurations at different times with small explicit symmetry breaking effects. The shading (color) representing 
dominant region in (d) corresponds to the true vacuum with = 0. (a) - (d) show the growth of domains for r = 1.2, 1.6, 2.0, 
and 2.4 fm (with corresponding values of temperature T = 376, 342, 317, 298 MeV). 



V. STRONG EXPLICIT SYMMETRY BREAKING AND LARGE FIELD OSCILLATIONS 

Now we consider the case when explicit symmetry breaking effects make the initial field configuration completely 
biased towards 9 — direction. Here the initial patch of I rolls down entirely towards 9 — direction with angular 
variations decreasing during the roll down. The dynamics of transition is entirely different in this case. Clearly there 
is no possibility of different Z(3) domains here, hence no Z(3) interfaces, or QGP strings will form. We also do not 
see any bubble like structures here as were seen in Fig. 2. Instead we find the I settles down to the true vacuum after 
undergoing huge oscillations, with large length and time scales. 

These large oscillations are very similar to coherent oscillations of the inflaton field in the context of inflation in 
the early Universe [7]. For the inflation, decay of such coherently oscillating field to particles can be via parametric 
resonance leading to novel features in the reheating of the Universe after inflation. Existence of similar oscillations 
here raises possibilities of parametric resonance for RHICE and similar novel dynamics of particle production during 
such early stages of QGP evolution. As for the universe, new possibilities of thermalization and symmetry changes 
may arise here. We hope to explore these issues in a future work. 

A direct implication of the existence of such huge oscillations of I in the context of RHICE is its possible effects 
on the growth of flow anisotropics. In the hydrodynamical studies of elliptic flow in non-central collisions it is known 
that much of the flow anisotropy develops during early stages of the plasma evolution 123] ■ In such studies 
one starts with equilibrium QGP phase where flow develops due to pressure gradients. In view of the possibility 
of large oscillations during early stages of transition to the QGP phase, the equilibrium starting point of these 
hydrodynamics simulations becomes suspicious. Our present study does not allow us to address this issue in the 
context of hydrodynamical evolution. However, even with the field theoretical simulation in this work, we can do a 
comparative study of momentum anisotropy development with and without the presence of large oscillations of the 
order parameter field I. For this we proceed as follows. 

First we need to model the initial QGP system with appropriate spatial anisotropies. One cannot then use the 
square lattice with uniform temperature. We use the temperature profile of Woods-Saxon shape with the size in the X 
and Y directions being different representing elliptical shape for a non-central collision. This allows us to have a well 
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defined size for the central QGP region, with temperature smoothly decreasing at the boundary of this region. The 
transverse size R of this system (i.e. profile of temperature) is taken to increase with uniform acceleration of 0.015 c 
per fm, starting from an initial value of R equal to the nuclear radius [53] . The initial transverse expansion velocity is 
taken to be zero. This expanding background of temperature profile is supposed to represent the hydrodynamically 
expanding quark-gluon plasma in which the evolution of the order parameter field / will be studied. It may appear 
confusing as I is expected to represent the QGP phase. Indeed, the normalization of the effective potential in Eq.(5) 
from refs.[5J 6 is carried out precisely so that it represents energy density and pressure of gluons plus quark degrees 
of freedom. Still, the dynamics of I from Eq.(5) does not carry the information of hydrodynamical degrees of freedom. 
Various particle modes in I need to be excited, which should reach equilibrium, and only then we can expect some 
sort of hydrodynamical evolution. Clearly the initial field configurations assumed here are far from representing such 
a hydrodynamical state. A consistent interpretation of our simulation can be that we are studying long wavelength 
modes of I which are couple to a background of short wavelength modes which are in thermal equilibrium. This 
equilibrated background is expanding with velocity as mentioned above, and drives the evolution of large wavelength 
modes of I via Eq(6). 

With this interpretation, our task is straightforward. The central temperature of the Woods-Saxon profile is taken 
to decrease by assuming that the total entropy (integrated in the transverse plane) decreases linearly as appropriate 
for Bjorken dynamics of longitudinal expansion. Note that, with the transverse expansion being non-zero now, the 
central temperature will decreases faster than r -1 / 3 . We show, in Fig. 4, a sequence of surface plots of the magnitude 
of I showing huge oscillations with large length scale during quench. Lattice here is again 2000 x 2000 but we take a 
large value of Ax = 0.25 fm so that the physical lattice size is 50 fm x 50 fm. The Woods-Saxon temperature profile 
(representing QGP region) is taken to have a diameter of about 16 fm as appropriate for Au-Au collision for RHICE. 
Large physical size of the lattice allows for the evolution of the QGP region to be free from boundary effects. 




FIG. 4: 

Surface plots of the magnitude of I with circular geometry for the QGP region, (a) - (d) show large oscillations of I during the 
quench when I rolls down everywhere towards the same vacuum with 8 = 0. Plots in (a)-(d) correspond to r = 2.6, 3.6, 4.4, 
6.2 fm/c, with corresponding values of temperature being T = 281, 251, 228, 191 MeV, respectively 

At any stage we can calculate the energy momentum tensor T^ v of I. We then calculate the spatial eccentricity e x 
of the I field configuration in the standard way, 



e = J dxdy(y 2 - x 2 )p 
x J dxdy(y 2 + x 2 )p 

Where p is the energy density. We can also calculate the momentum density at any time using T 0x and T 0y . Using 
these components we know the momentum density vector at every stage. By integrating it in angular bins we calculate 
its various Fourier coefficients, in particular the elliptic flow coefficient V2- 
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VI. EFFECTS OF LARGE I OSCILLATIONS ON FLOW ANISOTROPY 

To study the effects of large I oscillations on flow anisotropy, we consider two separate cases. Note that now we are 
considering the cases with explicit symmetry breaking, with its effect being strong on the initial field configuration. 
First we consider the quench case as described above. Here, we start with the initial field configuration corresponding 
to a small patch near the equilibrium point of T = effective potential. The patch is taken, as above, with field 
magnitude randomly varying in angle between and 2-7T with random amplitude uniformly varying from to 0.01 
times the vev of I but now shifted by a constant value of 0.011 times the vev along 6 — direction. This simulates 
the effect of strong explicit symmetry breaking as the entire patch rolls down towards 9 — direction. This patch 
is then evolved with the temperature profile of Woods-Saxon shape as described above with the central temperature 
having initial value equal to To = 400 MeV. As a sample case we give results for the case when the initial (elliptical 
shaped) temperature profile has an eccentricity of 0.5 (with the major and minor axes being along the x and y 
axes respectively), which for uniform energy density will correspond to e x — —0.143 from Eq.(7). The initial field 
configuration rolls down in the entire central region towards roughly = direction leading to strong I oscillations. 
e x , and flow coefficients, e.g. i>2 are calculated at each stage during the evolution. This is shown in Fig. 5. In this 
section we take the lattice to be 1000 x 1000 with the physical size of 25 fm x 25 fin. This still allows for sufficient 
separation of the QGP region (of size 16 fm diamater) from the boundary. Here we show plots (both for e x and for 
V2) for two different realizations of the random initial field configuration (shown as solid and dashed plots). Here, 
and in all the figures below, we will show plots for a time upto r = 10 fm/c, starting with r = tq = 1 fm/c at the 
initial stage. The central temperature decreases (now faster than r^ 1 / 3 due to nonzero transverse expansion) from an 
initial value T = T = 400 MeV at r = r = 1 fm/c to the final value T = 147 MeV at r = 10 fm/c. Note from the 
difference in solid and dashed plots in Fig. 5 that differences in the initial field configuration, which have very small 
magnitudes, lead to huge differences in the values of e x and i>2- 




-0.12 

-0.14 I 1 1 1 1 ' 1 1 1 

123456789 10 

x(fm) 
FIG. 5: 

These plots correspond to quench case when initial elliptical shaped temperature profile has eccentricity of 0.5. Solid and 
dashed plots correspond to two different realizations of the initial random field configuration, (a) and (b) show plots for e x 

and elliptic flow V2 respectively. 

This situation should be compared to the case of equilibrium initial conditions, as appropriate for conventional 
hydrodynamical simulations. This equilibrium initial condition is implemented in the following manner here. For 
the initial temperature profile of Woods-Saxon shape (again, with a given eccentricity), with initial temperature 
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T = Tb = 400 MeV, we first need to determine the appropriate I configuration which assumes vacuum expectation 
value everywhere depending on the local value of the temperature. To achieve this, we first evolve field configuration 
for 1000 time steps with highly dissipative dynamics while keeping the temperature profile to remain fixed at the initial 
profile. A reasonably smooth initial test profile of I is used as the initial field configuration. With highly dissipative 
evolution, the field everywhere settles to the local minimum of the potential quickly. The final configuration is found 
to be reasonably independent of the initial configuration assumed for I as long as it is smooth. This final configuration 
has correct profile as appropriate for the Woods-Saxon profile of temperature representing lowest energy configuration 
everywhere. In order to make a suitable comparison with the quench case we must incorporate small fluctuations 
around this equilibrium configuration everywhere. For this purpose we add to the value of the field everywhere, 
a small fluctuating field component with randomly varying angle between and 2-7T and with random amplitude 
uniformly varying from to 0.01 times the vev of I (as for the quench case). This new field configuration represents 
the equilibrium field configuration everywhere, with small fluctuations. This is taken as the initial field configuration 
for subsequent evolution where now any extra dissipation is switched off. This stage is taken as representing the initial 
time t = r = 1 fm. The field now evolves with the field equations, Eq.(6). Temperature profile also is now allowed 
to change in time as mentioned above. e x , and Vi etc. are calculated at each stage. Fig. 6 shows these plots for this 
equilibrium case, starting from the time slightly after when extra dissipation (to achieve equilibrium configuration) has 
been switched off. When random fluctuating component is introduced at every lattice site after the end of dissipative 
evolution, it introduces large gradients. Thus, for a very short time, there are large fluctuations as the field smoothens 
to certain level. Thus we show plots slightly after (by about 0.2 fm/c) the introduction of random field component. 




-0.06 1 — ' 1 ' 1 1 1 1 ' 1 

23456789 10 

T(fm) 



FIG. 6: 

Plots as in Fig. 5, now for the equilibrium case. Solid and dashed plots show different realizations of the small random 
fluctuating part of the equilibrium field configuration. The temperature profile has same initial eccentricity of 0.5 (with 
corresponding value of e x = —0.143) as in Fig.5. (a) and (b) show plots of e x and elliptic flow V2 respectively. Note, the initial 
value of e x for the equilibrated configuration in (a) is very close to -0.143 as shown by the plot in the inset. 

Comparison of Fig.5 and Fig. 6 shows the dramatic effects of quench induced oscillations on flow anisotropics. First 
note that for the equilibrium case the initial value of e x is close to the value -0.143 (as shown by the inset in Fig. 6) 
which exactly corresponds to the initial eccentricity of 0.5 for the temperature profile. This gives us confidence that 
our procedure of achieving equilibrated configuration works well. In contrast, initial value of e x for the case of quench 
in Fig.5 is very different showing the importance of fluctuations for this case. We further see huge fluctuations in the 
values of e x and i»2 in Fig.5 compared to the equilibrium case shown in Fig. 6. Note that though there are oscillations 
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in V2 for the equilibrium case also, they do not change much from one event to another. In contrast, for the quench 
case in Fig. 5, there are huge variations between the two different realizations, i.e. different events. These large 
fluctuations in Fig. 5 arise due to large oscillations in I which itself depends on the nature of randomness in the initial 
field configuration. Note that plots for e x in Fig. 6 show a sharp change near r = 5.5 fm after which e x settles down to 
a value close to zero. The central temperature at that stage is about 208 MeV. Presumably a large part of the region 
(away from the central part where temperature is lower due to a Woods-Saxon profile) may be undergoing transition 
at that stage leading to large changes in field dynamics. We have checked with the contour plots that the ellipticity 
of the profile of I, as well as that of energy density, does not undergo rapid changes during this stage. (In fact this is 
a reflection of a shortcoming in our model where the temperature profile is taken to have definite initial eccentricity, 
and its further expansion is with definite acceleration as discussed above, starting with zero transverse velocity. A 
more appropriate model may be to take the time evolution of the spatial eccentricity of the temperature profile from 
hydrodynamical models and study the evolution of I using that.) We note strong fluctuations in both these quantities 
around this stage, which may be responsible for these large changes in e x at this stage. We hope to develop better 
understanding of the dynamics during this stage and discuss in a followup work. At present what is important to 
note is that, apart from this peak region, everywhere else e x shows rather stable values settling down to a value close 
to zero after r ~ 5.5 fm, and does not fluctuate much, as compared to the quench case in Fig. 5a. 

To further illustrate the importance of fluctuations for the quench case, we show plots for the case with zero 
eccentricity of the QGP region (i.e. for the temperature profile) which also will mean zero value of e x for uniform 
density case. Thus, any non-zero e x and v-i arise only from the randomness in the initial field configuration. Fig. 7a 
and Fig. 7b show the plots for the equilibrium case. We note that e x and vi remain very small, apart from a large 
change near r = 5.5 fm for e x , just as in Fig. 6. Again, at present what is important to note is that, apart from this 
peak region, everywhere else e x and v% remain very small for equilibrium case, as expected for the zero eccentricity 
case. Different plots in Fig. 7a and Fig. 7b correspond to different realizations of the initial random field configuration. 
This situation should be contrasted with the quenched case as shown in Fig. 7c and Fig.7d. Note that e x and vi now 
fluctuate with large amplitudes, even though eccentricity of temperature profile is zero. Further, different realizations 
of the initial random field configurations lead to widely different plots of these quantities. This shows the fluctuating 
nature of development of flow anisotropics if large I fluctuations are present. These results suggest the dynamics of 
the order parameter field, especially such large length scale oscillations may play an important role in determining 
flow anisotropies and it needs to be incorporated in the hydrodynamical models. 

VII. CONCLUSIONS 

We carry out 2+1 dimensional simulation to study the dynamics of confinement-deconfinement transition as a 
quench for the central rapidity region in relativistic heavy-ion collision experiments. We work in the framework of the 
Polyakov loop model. The initial field (in the confining phase) is taken to cover a small neighborhood of the confining 
vacuum I ~ as appropriate for the initial T = system. This initial field I, unable to relax to the new equilibrium 
vacuum state with quenched potential at T = T = 400 MeV, becomes unstable and rolls down in different directions 
from the top of the central hill in the effective potential of I. We study the formation of Z(3) domain structure during 
this evolution. When explicit Z(3) symmetry breaking effects (arising from dynamical quark effects) are small, then 
we find well defined Z(3) domains which coarsen in time. During early stages Z(3) domains (and Z(3) domain walls) 
have small sizes of order 1 fm. However at this stage the magnitude of I is very small, of order few present of its 
vacuum expectation value. Domains coarsen rapidly as the magnitude of I grows and by the time r = 4 — 5 fm, 
domains are of size several fm. Surprisingly, the magnitude plot of I shows vacuum bubble like configurations, such 
as those which arise in a first order transition. This first order transition like behavior occurs even though there is 
no metastable vacuum separated by a barrier from the true vacuum for the parameter values used. These bubble like 
configurations expand as well, somewhat in a similar manner as during a first order transition. This result points to 
new interesting possibilities for the phase transition dynamics. For example, such bubble like structures may lead to a 
dynamics of phase separation in the case of the universe similar to the original Witten's scenario [20j . even when there 
is no underlying first order transition. This may also have important implications for RHICE. This new possibility of 
transition dynamics needs to be understood and analyzed in detail, see [5T] for a study of these issues. 

When the initial patch of I is only partially symmetric around I = (as appropriate for small explicit symmetry 
breaking from quark effects), the dynamics retains these qualitative aspects, with true vacuum domains (with 8 = 0) 
growing dominantly at the cost of the other two metastable Z(3) domains. In this case Z(3) walls are fewer and 
relatively smaller, and they disappear more quickly. When the initial patch of I (around equilibrium point for T = 
effective potential) is significantly shifted towards the true vacuum for the quenched T = T effective potential (as will 
happen when explicit symmetry breaking is strong), then I rolls down roughly along the same direction with angular 
variations becoming smaller during the roll down. In this case only 9 = vacuum survives and no other Z(3) domains 
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FIG. 7: 

These plots are for the zero eccentricity for the temperature profile. In all the figures here, solid and dashed plots correspond 
to different realizations of initial random field configuration, (a) and (b) show plots of e x and i>2 for the equilibrium transition 

case, while (c) and (d) correspond to the quench case. 
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are formed. Also, in this case we do not find bubble-like configurations. However, in this case we find huge oscillations 
of I with large length scales. This is similar to the scenarios of reheating via parametric resonance in the case of 
inflation in the early Universe. In our case of RHICE also it raises important questions about the possibility of novel 
modes of particle production from these large oscillations of I during the early stages of the transition. We have shown 
that these large I oscillations can strongly affect the evolution of flow anisotropics in RHICE. The spatial eccentricity 
and the flow coefficients, e.g. the elliptic flow V2 are found to undergo large fluctuations during the evolution of the 
system. These results suggest the dynamics of the order parameter field, especially such large length scale oscillations 
may play an important role in determining flow anisotropies and it needs to be incorporated in the hydrodynamical 
models. 
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